PL— TR— 94— 3037 


PL-TR- 

94-3037 


AD-A283  968 

ill 


Engineering  Considerations  for  the  Self- Energizing  Magneto  - 
plasmadynamic  (MPD)  -  Type  Fusion  Plasma  Thruster 


Chan  K.  Choi 
Larry  T.  Cox 


Purdue  University 
School  of  Nuclear  Engineering 
\1290  Nuclear  Engineering  Building 
West  Lafayette  IN  47907—1290 


August  1994 


Final  Report 


^ 94-27941 

IliililHIl 


g  29  2  35 

PHILLIPS  LABORATORY 
Propulsion  Directorate 

AIR  FORCE  MATERIEL  COMMAND 

EDWARDS  AIR  FORCE  BASE  CA  93524-7001 


NOTICE 


When  U.S.  Government  drawings,  specifications,  or  other  data  are  used  for  any  purpose 
other  than  a  definitely  related  Government  procurement  operation,  the  fact  that  the 
Government  may  have  formulated,  furnished,  or  in  any  way  supplied  the  said  drawings, 
specifications,  or  other  data,  is  not  to  be  regarded  by  implication  or  otherwise,  or  in  any  way 
licensing  the  holder  or  any  other  person  or  corporation,  or  conveying  any  rights  or  permission 
to  manufacture,  use  or  sell  any  patented  invention  that  may  be  related  thereto. 

FOREWORD 

This  report  was  prepared  by  Purdue  University,  West  Lafayette  IN,  under  contract 
F04611-90-K-0054,  for  Operating  Location  AC,  Phillips  Laboratory,  Edwards  AFB,  CA 
93524—  7048.  Project  Manager  for  Phillips  Laboratory  was  Dr.  Franklin  B.  Mead. 

This  report  has  been  reviewed  and  is  approved  for  release  and  distribution  in  accor¬ 
dance  with  the  distribution  statement  on  the  cover  and  on  the  SF  Form  298. 


Project  Manager 


Chief,  Emerging  Technologies  Branch 


Acting  Director, 

Fundamental  Technologies  Division 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No  0704-0188 

Pvbttc  rtportia|  burdca  for  Ihb  collection  o f  brformatioa  Is  estireaUd  to  average  1  boor  per  response,  including  the  time  for  reviewing  instructions 
searching  existing  data  sources  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send 
comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information.  Including  suggestions  for  reducing  this  burden  to 
Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington, 

VA  22202-4302,  and  to  the  Office  of  Maaacemeat  and  Budget  Paperwork  Reduction  Protect  (0740-0188),  Washington  DC  20503, 

1.  AGENCY  USE  ONLY  (LEAVfe  BLANK  2.  REPORT  DATE  3,  REPORT  TYPE  AND  DATES  COVERED 

August  1994  Final  July  1992  -  Dec  1993 

1  ir[gfn^nngBPonlideniUons  for  the  Self-Energizing 
Magnetoplasmadynamic  (MPD)  -  Type  Fusion  Plasma 

Thruster 

5.  FUNDING  NUMBERS 

c,  F0461 1-90-K-0054 

PE.-  62302F 

PR:  3058 

TA:  00 AF 

6.  aUTH6R(S)  ““ 

Chan  K.  Choi 

Larry  T.  Cox,  Jr. 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Purdue  University 

School  of  Nuclear  Engineering 

1290  Nuclear  Engineering  Building 

West  Lafayette  IN  47907-1290 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Phillips  Laboratory 

OL-AC  PL/RKFE 

9  An  tares  Road 

Edwards  ABF  CA  93524-7860 

10.  SPONSORING/MONITORING 

AGENCY  REPORT  NUMBER 

PL-TR-94-3037 

11.  SUPPLEMENTARY  NOTES 

COSATI  CODE(S):  2013 

12a.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  Public  Release;  Distribution  is  Unlimited 

12b.  DISTRIBUTION  CODE 

A 

13.  ABSTRACT  (MAXIMUM  200  WORDS) 

With  a  Dense  Plasma  Focus  (DPF)  device  as  the  central  concept,  studies  have  been  done  to  determine  its 
feasibility  as  a  propulsion  system  for  space.  In  this  report,  past  work  in  the  areas  of  propulsion  system  code  development 
is  discussed,  as  well  as  the  recent  work  dealing  wit  stability  analysis  and  scaling  laws  in  the  pinch  region.  A  modeling 
based  on  a  tokamak-like  m=0  instability  relating  to  the  electron  drift  velocity  wavelength  is  established.  Magnetic  Held 
and  kinetic  temperature  profiles  are  calculated  based  in  certain  assumptions  about  the  plasma.  It  is  found  that  the 
results  of  the  pinch  equilibrium  profiles  agree  with  the  assumed  current-squared  scaling  of  kinetic  temperature  density, 
and  that  the  total  fusion  power  released  from  the  pinch  scales  as:  Pf  ~  Ip5*4*.  This  is  found  to  correlate  well  with  the 
reaction  rate  parameter  data  for  the  utilized  ftiel  of  deuterium  and  helium-3.  The  resulting  profile  shape  for  kinetic 
temperature  agrees  in  form  with  the  proposed  modeling.  Using  the  more  detailed  pinch  calculations  in  the  already 
existent  DPF  propulsion  code,  it  is  found  that  for  optimal  performance  a  current  range  of  30  to  40  [MA]  is  needed  to 
obtain  results  found  in  earlier  work. 

14.  SUBJECT  TERMS 

fusion  propulsion;  dense  plasma  focus;  magnetoplasmadynamic  thruster;  MPD 
thruster;  advanced  fuel  (D-3He)  fusion 

15.  NUMBS*  OF  PAGES 

It.  PRICE  CODS 

i7.  security  aAMncAnoN  it.  sicimrrr  classification  <’•  sscvsrrv  classification 

orsarosT  ofthispaci  op  abstract 

Unclassified  Unclassified  Unclassified 

2t  LIMITATION  OP  ABSTRACT 

SAR 

i/ii 


NSN  7540-010280-5500 


Standard  Form  298  (Rev  2-89) 
Prescribed  by  ANSI  Std  239-18 
298-102 


ABSTRACT 


With  a  Dense  Plasma  Focus  (DPF)  device  as  the  central  concept,  studies  have  been  done  to 
determine  its  feasibility  as  a  propulsion  system  for  space.  In  this  report,  past  work  in  the  areas  of 
propulsion  system  code  development  is  discussed,  as  well  as  the  recent  work  dealing  wit  stability 
analysis  and  scaling  laws  in  the  pinch  region.  A  modeling  based  on  a  tokamak-like  m=0  instability 
relating  to  the  electron  drift  velocity  wavelength  is  established.  Magnetic  field  and  kinetic 
temperature  profiles  are  c. '  T»H  based  in  certain  assumptions  about  the  plasma.  It  is  found  that 
the  results  of  the  pinch  es  1  «  profiles  agree  with  the  assumed  current-squared  scaling  of 

kinetic  temperature  density,  and  ‘hat  the  total  fusion  power  released  from  the  pinch  scales  as:  Pf  ~ 
Ip5-4-.  This  is  found  to  correlate  .veil  with  the  reaction  rate  parameter  data  for  the  utilized  fuel 
of  deuterium  and  helium-3.  The  resuu»ng  profile  shape  for  kinetic  temperature  agrees  in  form 
with  the  proposed  modeling.  Using  the  more  detailed  pinch  calculations  in  the  already  existent 
DPF  propulsion  code,  it  is  found  that  for  optimal  performance  a  current  range  of  30  to  40  [MA] 
is  needed  to  obtain  results  found  in  earlier  work. 
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INTRODUCTION 


In  this  third,  final  report  of  the  three-year  project,  primary  work  is  centered  on  the 
equilibrium  and  stability  analyses  of  pinch  plasmas  and  the  scalings  of  dense  plasma  focus  (DPF) 
properties  with  respect  to  various  parameters.  The  primary  properties  of  interest  in  the  pinch 
scalings  are  the  plasma  pinch  length  and  radius  and  the  fusion  product  yield.  A  model  based  on  a 
tokamak-like  m=0  sausage  instability  relating  to  the  electron  drift  velocity  wavelength  was 
developed  as  described  in  the  later  section.  Regarding  the  length  and  radius  of  the  pinch  region, 
the  parameters  of  interest  have  been  the  average  ion  current,  the  plasma  kinetic  temperature  within 
the  pinch,  and  the  plasma  number  densities  within  the  pinch.  For  the  fusion  product  yield,  it  has 
been  investigated  how  it  varies  with  the  ion  current  flowing  in  the  pinch  region  of  the  plasma  focus 
structure. 

For  the  properties  of  pinch  length  and  radius,  a  model  was  developed  based  on  the  electron 
drift  velocity.  At  its  heart  is  the  coupling  of  the  wave  number  of  the  electron  drift  velocity  and  the 
wave  number  associated  with  the  ion  gyromagnetic  drift  (the  Larmor  radii).  A  variant  form  of  the 
traditional  cylindrical  coordinate  system  used  with  the  DPF  has  been  developed  to  illustrate  how 
this  model  is  to  be  applied. 

With  regard  to  the  fusion  product  yield,  a  two-dimensional  grid  was  used  to  divide  the 
pinch  region  into  divisions  along  the  traditional  r  and  z  coordinates.  Assuming  an  axial  velocity 
profile  and  an  initial  number  density  profile  along  these  axes,  as  well  as  a  parabolic  contour  for  the 
edge  of  the  pinch  region,  the  ion  kinetic  temperature  profiles  were  found.  This  was  done  based  on 
the  magnetic  field  penetrating  to  the  centerline  of  the  pinch  before  going  to  zero.  It  was  assumed 
that  the  total  pressure  from  the  sum  of  the  magnetic  field  pressure  and  the  plasma  particle  pressure 
remained  constant  throughout  the  region  of  the  pinch.  After  determining  the  kinetic  temperature 
profiles,  the  fusion  product  yield  is  found  based  on  the  reaction  rate  parameter  data  and  the 
calculated  volume  of  the  pinch  region.  This  may  be  correlated  with  the  calculated  current  within 
the  pinch  region  in  order  to  determine  the  scaling  of  the  fusion  product  yield  with  varying  pinch 
current. 

Use  was  made  of  the  codes  developed  earlier  in  project  history  [1,2]  in  order  to  determine 
the  effects  of  the  recent  work  on  propulsion  performance  parameters  for  various  mission 
specifications.  The  effects  on  mission  parameters  such  as  velocity  increment  (Av),  propellant  mass 
flow  rate,  number  of  thrusters  used,  thrust-to-weight  (F/W)  ratio,  and  specific  impulse  (Isp)  have 
been  investigated. 
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PRELIMINARY  ANALYSIS  AND  REVIEW  OF  PREVIOUS  WORKS 


The  dense  plasma  focus  (DPF)  is  one  of  a  class  of  magnetically-confined  fusion  (MCF) 
devices  known  as  Z-pinches.  One-dimensional  (1-D)  Z-pinch  devices  are  characterized  by  a 
"straightened"  toroidal  geometry  with  a  purely  poloidal  magnetic  field.  They  possess  a  toroidal 
geometry  from  the  standpoint  that  the  plasma  behavior  is  quite  similar  to  that  which  would  be 
expected  in  a  single,  extracted  section  of  a  torus  device.  Using  cylindrical  coordinates,  as  shall  be 
the  system  of  choice  throughout  this  report,  only  the  Be  component  of  the  magnetic  field  B  will  be 
considered  •>  be  non-zero.  It  is  induced  by  a  longitudinal  current  Iz  flowing  in  the  plasma.  When 
the  plasma  current  vanishes,  no  background  B  field  remains.  The  following  arguments  show  the 
derivation  of  the  radial  pressure  balance  in  a  1-D  Z-pinch  [3]. 

Assuming  azimuthal  symmetry  of  B,  the  equation: 


VB  =0 


(2.1) 


appears  as: 


1  aBe 

r  30 


=  0  . 


(2.2) 


Applying  Ampere's  Law  shows  that  Jz(r)  is  the  only  non- zero  component  of  the  current  density: 


J.-m.'tt (rBe>  • 


Also,  only  the  radial  component  of  the  momentum  equation  is  nontrivial: 

JB 


(2.3) 


dr 


(2.4) 


If  Eq.  (2.3)  is  substituted  into  Eq.  (2.4): 


dp 

dr 


+ 


Be  d 

M^dr 


(rB„)  =  0  , 


(2.5) 
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which  may  also  be  written  in  the  form: 


df 


+ 


2  \ 


Be_ 

2U0 


B 

40r 


=  0  . 


(2.6) 


Eq.  (2.6)  is  the  basic  radial  pressure-balance  relation  for  a  1-D  Z-pinch.  The  two  terms  within  the 
parentheses  represent  the  particle  pressure  and  the  magnetic  pressure,  respectively.  The  final  term 
represents  the  tension  force  generated  by  the  B  field  line  curvature.  This  tension  force  is  of  utmost 
importance  in  a  Z-pinch  device,  as  its  presence  provides  the  radial  confinement  of  the  plasma. 
Later,  the  DPF  will  be  studied  from  a  two-dimensional  standpoint  in  the  r  and  z  coordinates, 
moving  from  the  1-D  argument  in  r  presented  above.  Azimuthal  symmetry  will  continue  to  be 
assumed  to  account  for  the  0  component. 

The  Dense  Plasma  Focus  (DPF)  was  first  engineered  in  separate  endeavors  by  Filippov  [4] 
in  1962  and  Mather  [5]  in  1964.  Most  of  the  DPF  experimental  devices  fall  into  one  of  two 
primary  geometries,  accompanied  by  the  corresponding  name  of  one  of  these  two  initial 
researchers.  The  Mather-type  plasma  focus  device  uses  axial  injection  and  radial  compression  of 
the  plasma  to  form  the  pinch,  employing  an  electrode  geometry  having  an  aspect  ratio  a  <  1,  while 
the  Filippov  device  is  characterized  by  a  >  1.  The  aspect  ratio  is  defined  as  the  ratio  of  cathode 
diameter  to  anode  length.  A  schematic  comparison  is  shown  in  Figure  2.1.  The  DPF  combines 
aspects  of  two  other  Z-pinch  designs:  the  compressional  Z-pinch  and  the  gas-embedded 
(imploding  liner)  Z-pinch.  It  forms  in  much  the  same  manner  as  the  compressional  pinch,  but  has 
higher  plasma  densities  on  the  order  of  those  found  in  the  gas-embedded  pinch  [6].  More  than  40 
laboratories  in  Europe,  the  United  States,  Russia,  India,  and  Japan,  among  other  nations,  currently 
perform  DPF  experimental  and  theoretical  research  [7]. 

Of  utmost  interest  to  the  success  of  a  DPF  device  as  a  power  source  and  a  propulsion 
device  is  the  production  of  energetic  particles  from  the  thermonuclear  processes  which  occur  in  the 
pinch  region  of  the  compressed  plasma.  Most  experimental  work  done  to  date  has  used  a  fusion 
fuel  mixture  of  deuterium.  Deuterium  has  the  advantages  of  being  both  plentiful  (and  therefore 
inexpensive)  and  non-radioactive,  as  opposed  to  the  usage  of  tritium  in  conjunction  with 
deuterium. 

The  deuteron  (D)  -  deuteron  (D)  reaction  occurs  through  two  primary  pathways: 

D+  +  D+  — >  p+  +  T+  +  4.03  MeV 
D+  +  D+  — >  n  +  3He++  +  3.27  MeV 
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Figure  2.1 

Mather  vs.  Fil-ppov  Plasma  Foci  Aspect  Ratios  [7] 


Because  the  DPF  has  long  been  used  as  a  laboratory  neutron  source,  one  of  the  primary  goals  of 
research  has  been  to  maximize  neutron  yield.  This  necessarily  indicates  a  large  thermonuclear 
reaction  rate.  In  working  to  maximize  neutron  production,  it  was  necessary  to  determine  which 
parameters  of  DPF  operation  affected  the  neutron  yield  and  to  what  degree.  One  of  the  key 
parameters  in  this  work  has  been  the  current  within  the  plasma  pinch.  Much  has  been  done  from 
both  a  theoretical  and  an  experimental  standpoint  to  determine  the  degree  of  effect  these  currents 
have  on  the  neutron  yield. 

Some  of  the  first  theoretical  scaling  work  was  done  by  Kaeppeler  in  1974  and  1976  [8,9], 
who  determined  that  the  neutron  yield  Yn  should  scale  with  the  pinch  plasma  current  as  Ip4. 
Decker  deduced  in  1977  [  10]  a  scaling  relation  of  Yn  ~  Imaxb,  with  Imax  being  the  maximum  main 
bank  capacitor  discharge  current  and  the  exponent  b  varying  between  3.5  and  4.6  .  Table  I  lists  a 
summary  of  the  theoretical  work  found  in  the  literature. 


Table  I 

Researchers  with  Theoretical  Neutron  Yield  Scaling  Proposed 


Researcher/Year 

LocjitiQn/Dcvice 

Yu 

Kaeppeler/ 1974[8] 

Stuttgart 

V 

Rapp/1974[12] 

Stuttgart/Nessi 

I  4 

AP 

Decker/ 1977[  11] 

Stuttgart 

T  3.5  to  4.6 

Amax 

Moving  to  the  experimental  front,  in  1978  Shyam  and  Srinivasan  found  using  a  1(H)- Joule 
DPF  that  Y„  scaled  with  W1 73  and  I4  29  [13],  where  W  is  the  capacitor  bank  energy  and  I  is  the 
discharge  current  of  the  bank.  Mather  found  Yn  -  I4  in  1971  [14]  and  Bernard  found  Yn~  I3  3  in 
1976  [15,16],  based  on  their  experimental  findings.  In  more  recent  years,  scaling  law 
development  has  focused  on  the  use  of  the  plasma  current  Ip  in  scaling  relations,  rather  than  the 
discharge  current  I.  This  is  due  to  the  fact  that  the  current  flowing  within  the  pinch  is  directly 
characteristic  of  the  pinched  region,  where  the  fusion  reactions  occur,  and  should,  therefore,  be 
closely  interwoven  with  the  mechanisms  involved  with  fusion  particle  production  in  the  pinch. 
Scaling  with  the  discharge  current  fails  to  account  for  the  current  loss  during  the  rundown  phase  of 
the  plasma  sheath  which  forms  the  pinch.  Recognizing  that  the  presence  of  such  current  losses 
leaves  a  degree  of  uncertainty  in  the  number  of  particles  trapped  in  the  pinch  and  that  improvements 
in  measuring  pinch  plasma  currents  in  recent  years  have  taken  place,  the  scaling  laws  have  surely 
been  directed  away  from  using  the  discharge  current  as  a  key  parameter  when  considering  fusion 
yield.  Experiments  have  shown  Y„  ~  Ip4  1  for  the  Minifokus  device  at  Stuttgart  in  1977  [  17]  and 
Conrads  has  found  Y„  ~  Ip3  3  based  on  agreement  of  the  SPEED  II  device  at  the  Institut  fur 
Plasmaphysik  in  Germany  with  worldwide  research  in  recent  years  [6],  Stygar,  et  al.  determined  a 
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scaling  of  Ip4  4  in  1982  [18].  A  summary  of  these  experimentally  determined  scalings  are 
presented  in  Table  II,  along  with  some  device  characteristics  of  those  used  by  the  researchers 
indicated  below. 

Of  other  recent  interest  has  been  the  study  of  using  the  combination  of  the  fusion  fuels, 
deuterium  and  helium-3  [1,2,23-27],  as  opposed  to  the  traditional  D-D  combination  used  in 
experimental  facilities  and  the  radioactive  D-T  combination  used  in  tokamak  facilities  in  the  recent 
past.  The  advantages  of  such  a  fuel  combination  are  illustrated  by  observing  the  fusion  products 
produced  in  the  reaction: 


D+  +  3He++  — >  p+  +  4He++  +  18.3  MeV, 

where  the  products  are  evidently  only  charged  particles.  Such  reaction  products  have  the 
advantage  of  being  directed  using  a  magnetic  nozzle,  a  characteristic  which  lends  itself  for  use  in  a 
propulsion  system,  where  uni-directional  thrust  is  desirable.  Also,  looking  to  future  technological 
operations,  the  presence  of  charged  ions  makes  such  a  fuel  combination  attractive  for  use  in  a  direct 
conversion  scheme  to  produce  electricity  [28].  The  possible  D-D  background  fusion  reactions  are 
considered  to  be  suppressible  by  polarizing  the  spin  of  the  interacting  deuterons  [29]. 

Turning  to  the  scaling  of  the  pinch  size  (radius  and  length)  itself,  some  of  the  most  current 
work  has  been  done  by  Gerdin,  et  al.  using  a  device  at  the  University  of  Illinois  [30],  During  the 
late  1980's  they  measured  the  pinch  current  Ip  and  filling  pressure  p„  in  order  to  determine  their 
scaling  relation  with  the  plasma  pinch  length  h,  in  keeping  with  the  original  notation.  The  results 
of  the  pressure  scaling  measurements  are  included  for  the  sake  of  completeness,  but  the 
investigation  of  the  pressure  dependence  is  not  a  primary  focus  of  this  report.  Resulting  from  the 
work  of  Gerdin,  etal.  is  the  relation: 


•  (2.7) 

where  x  =  -0.02  ±0.12  and  y  =  0.26  ±0.13.  The  dependence  of  the  plasma  radius  r0  was  found 
to  relate  as: 


r„  - 1;  pj  .  (2.8) 

where  a  =  -0. 18  ±  0. 1 1  and  b  =  0.2  ±  0.03.  The  graphs  showing  plasma  length  and  radius  versus 
the  plasma  current  will  be  reproduced  in  a  later  section  on  the  pinch  scaling  relations.  Also  present 
in  Ref.  [30]  is  a  comparison  of  other  devices  where  the  value  of  h  has  been  predicted  using  the 
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scaling  law  presented  in  Eq.  (3.1).  The  pertinent  device  locations  studied,  along  with  the  results, 
are  found  in  Table  III. 


Table  II 

Researchers  with  Device  Parameters  and  Experimentally-Determined 
Neutron  Yield  Scaling  Relations  (NS=None  Stated) 


Researchers) 

Aulhfllfsl 

Bemard[15,16] 

Device 

(Location) 

(Limeil) 
Focus  1 
Focus  2 
Focus  3 
Focus  4 
Focus  5 
Focus  6 

wiM) 

15 

16 

30 

32 

32 

48 

Po 

(ton) 

Decker[ll] 

(Conrads[6]) 

SPEED  II 
(Stuttgart) 

1000 

0.1-10 

Decker[  19,20] 

Nessi 

(Stuttgart) 

30-60 

Decker[ll] 

HV-Focus 

12 

8 

Herold[21] 

Poseidon 

(Stuttgart) 

280-500 

4-6 

PF-360 

(Stuttgart) 

200 

12 

Mather[14] 


Schmidt[7]  Various 


Shyam/ 

TPF-1 

0.41- 

Srinivasan[22] 

(Bombay) 

2.15 

Shyam/ 

TPF-1 

0.1  0.1-10 

Srinivasan[13] 

(Bombay) 

Stygar[18] 

(Illinois) 

6-12.5  3 

Oppenlander 

Minifokus  12 

1.9 

[17] 

(Stuttgart) 

lAlcm}  1£(CT)  rptenr  Yqjt 

P  3 


I  3.3 

*P 

I4 


1.20 

3.0 

NS 

6.55- 

10.4 

5.0- 

6.0 

13.3- 

14.6 

7.5- 

11.25 

T  3.3 

Amax 

T  3.2 
1  m  ax 

I4 

13.3 

I  4 

Ap 

0.2- 

0.26 

0.8- 

1.8 

NS 

1.1 

3.6 

14.29 

J4.4±0.3 
(Peak  Y) 
j4.6±0.3 

I4  6(hot 
target) 
I52(cold 
target) 

2.5  4.5  I„4-i 
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Table  III 

Comparison  of  Plasma  Pinch  Length  h  in  Various  Mather-type  DPF  Experiments 
Using  Empirical  Scaling  Relation  (Eq.(3. 1))  [30] 


Location  and/or  Device 

rA  (cm) 

IfilkAl 

Da  (torr) 

fn  (cm) 

Exper. 
h  (cm) 

Eq.(3.1) 
h  (cm) 

U.  of  Illinois 

2.5 

550 

3.0 

0.15 

1.3 

1.3 

Stuttgart  Nessi 

3.3 

290 

1.5 

0.15 

2 

1.3 

Swierk 

5.0 

1000 

1.5 

0.38 

2-2.5 

3.3 

Darmstadt 

0.8 

200 

2.6 

0.04 

0.25 

0.35 

Stevens  Inst.  Tech. 

1.7 

400-500 

8 

0.2 

1.0  -  1.5 

1.7 
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EQUILIBRIUM  PROFILES  AND  STABILITY  ANALYSIS 


By  studying  the  equilibrium  profiles  of  the  pinch  properties,  one  can  determine  the  effects 
of  pinch  current,  pinch  radius,  pinch  length  and  other  parameters  on  such  properties. 
Determination  of  the  relations  between  such  properties  as  pressure,  plasma  temperature,  and  local 
magnetic  field  strength  with  pinch  radius  and  length  allows  one  to  see  how  such  properties  can  be 
motivated  to  produce  a  pinch  which  yields  the  maximum  thermonuclear  yield  of  particles  and 
energy. 

As  an  attempt  at  relating  the  pinch  radius  and  length  to  another  property,  one  may  consider 
modeling  the  pinch  length  as  one-half  of  the  wavelength  of  the  m  =  0  instability  which  forms  the 
pinch.  This  is  illustrated  in  Figure  3.1.  The  m  =  0  instability  presents  a  bridge  of  similarity  to  a 
tokamak  device,  where  disruptions  can  occur  along  the  toroidal  plasma  due  to  the  uneven  magnetic 
field  profile.  Now,  to  assess  the  feasibility  of  the  aforementioned  model,  equilibrium  profiles  need 
to  be  assumed.  If  such  profiles  are  assumed  to  behave  in  a  similar  manner  to  those  of  a  tokamak 
device,  which  is  also  a  device  of  toroidal  geometry  exhibiting  m  =  0  instabilities,  then  the 
following  relations  are  plausible: 


n(r)  _  kTe(r) 
nG  LTeo 


kTi(r) 

kTi,0 


=  1  -  0.9 


r 

r2 

‘max 


(3.1) 


where  n  is  the  number  density,  kT  is  the  kinetic  temperature  of  the  plasma,  and  r  is  the  radius  of 
the  pinch,  with  rmax  being  the  maximum  pinch  radius.  For  this  calculation,  the  coordinate  system 
shown  in  Figure  3.2  is  used.  In  order  to  have  the  same  arrangement  of  magnetic  field  as  in  a 
tokamak,  the  DPF  region  is  viewed  as  a  cross  section  of  the  tokamak  torus.  The  traditional  r-axis 
of  the  pinch  is  relabeled  the  z-axis,  which  is  shown  coming  out  of  the  page.  The  magnetic  field 
also  is  coming  out  of  the  page  at  the  upper  edge  of  the  pinch  region.  The  traditional  direction  of  the 
0-axis  has  been  rotated  90  degrees,  as  shown  in  Figure  3.2,  which  the  r-axis  used  in  Figure  3.2 
represents  the  distance  from  the  center  of  the  pinch  to  any  given  point  within  or  at  the  edge  of  the 
pinch.  The  curves  extending  behind  and  to  the  left  of  the  pinch  represent  the  tokamak  torus  which 
is  used  as  the  basis  for  coupling  the  drift  instability  to  the  DPF  pinch.  It  is  further  assumed  that 
rinax  is  constant  for  a  given  z  in  the  new  coordinate  system,  although  in  reality  it  would  be  a 
function  of  6  in  the  new  coordinate  system. 

Now,  since  particle  pressure  is  given  by  p  =  nkT,  the  electron  pressure  is  given  by: 
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Pinch  Length  =  1/2  of  m=0  Instability  Wavelength 


Modeling  of  Plasma  Pinch  as 


Figure  3.1 

One-Half  Electron  Drift  Velocity  m=0  Instability  Wavelength 
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Figure  3.2 

Hybrid  Coordinate  System  for  Purposes  of  Drift  Instability  Calculations 


(3.2) 


( 

Pe  nkTe  —  noekTt)e 

v 


r2  f 

1  -  0.9— 

^max  / 


According  to  Dolan  [31],  the  electron  drift  velocity  is  given  by: 


(3.3) 


in  terms  of  the  present  hybrid  coordinate  system.  The  wave  number  associated  with  this  drill 
velocity  matching  the  wave  number  associated  with  the  ion  cyclotron  frequency  is  given  by: 


(3.4) 


where  COci  is  the  cyclotron  frequency  and  A*  is  the  wavelength  of  the  drift  wave.  Combining  Eqs. 
(3.2-3.4)  yields  the  expression  for  A*.  using  0)ci  =  eB/mi: 


7.2711X1^7,3,. 

e2W 


(3.5) 


Next,  it  is  assumed  that  the  sum  of  the  particle  kinetic  pressure  and  the  confining  magnetic 
pressure  is  constant: 


B  l  B2v 

p  +  _.conSi  =  _  ,  0.6) 

where  p0  is  the  permeability  of  free  space,  4tcT0'7  [H/m],  and  Bov  is  the  value  of  the  magnetic 
field  at  the  vacuum  boundary.  Now,  Eq.  (3.6)  is  solved  for  Bz2  and  coupled  with  Eq.  (3.5), 
yielding: 


2.93’IQ19 

^max^oe 


(3.7) 


where  the  mj  value  used  is  the  mass  of  a  deuteron  in  order  to  compare  with  experimental  devices, 
which  use  deuteron-deuteron  fusion  fuel  mixtures.  Since  the  length  of  the  pinch  (/p)  is  equal  to 
one  half  of  and  assuming  values  of  lp  =  0.013  [m]  and  rmax  =  0.0015  [m]  [30],  the  value  of  noe 
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is  found  to  be  7.54- 1023  [nr3].  This  is  -  l()24  [nr3],  which  is  on  the  order  of  those  number 
densities  found  in  Z-pinches  by  Choi  of  2- 1024  [nr3]  and  in  dense  plasma  foci  by  Herold,  ei  til.  of 
3- 1024  [nr3]  [7].  Increases  of  one  to  two  orders  of  magnitude  in  number  density  will  be  required 
for  space  propulsion  applications,  as  will  be  made  evident  by  observing  the  values  of  parameters  in 
test  cases  shown  later.  It  is  generally  accepted,  however,  that  present  devices  can  readily  attain  ion 
number  densities  of  1025  [nr3]. 

Building  on  this  novel  approach  to  pinch  dimensions,  attention  is  turned  to  the  range  of 
influence  of  the  magnetized  plasma  in  the  pinch.  The  drift  is  found  to  proceed  along  the 
approximately  elliptical  pinch  boundaries  in  a  rough,  counter-clockwise  motion.  This  drift 
motion,  in  turn,  induces  a  magnetic  field  which  opposes  the  applied  Be  field.  This  is  known  as  the 
diamagnetic  effect.  The  diamagnetism  decreases  the  magnetic  field  as  r  approaches  the  central 
pinch  axis.  The  qualitative  significance  of  this  will  be  illustrated  later.  To  see  the  degree  of  effect 
diamagnetism  has  in  the  DPF,  the  Larmor  radius  of  the  deuterons  in  the  plasma  is  calculated,  as 
well  as  the  skin  depth.  The  Larmor  radius  frO  is  given  by: 


Tl~  ZeB 


(3.8) 


where  'i'  denotes  the  form  for  ions.  Here,  the  mass  mi  is  the  average  mass  of  a  deuteron  and  a 
helium-3  ion,  while  kT,  is  the  kinetic  temperature,  or  energy,  of  both  ion  species.  A  plasma  of 
these  two  components  requires  an  ignition  temperature  of  -  30  [keV],  so  this  value  is  chosen  for 
the  calculation  [24,25].  For  a  deuteron  the  atomic  number  Z  =  1,  while  Z  =  2  for  a  helium-3  ion, 
so  Z  =  1.5  is  used.  The  electronic  charge  e  is  1.602- 1019  [C].  To  determine  the  value  of  B 
needed,  a  balance  of  the  magnetic  compressional  pressure  and  the  particle  inertial  pressure  is 
assumed,  i.  e.,  the  commonly  referred  to  "(J  =  1"  condition,  where: 


P  = 


Ppart 

Pmag 


(3.9) 


and  the  summation  is  over  the  s  many  species  of  ions  plus  the  electrons  as  another  species,  with 
Ppart  representing  the  particle  pressure  and  pmag  representing  the  confining  magnetic  pressure. 
Using  the  number  density  n  =  7.54- 1023  [m*3]  found  using  Eq.  (3.7),  the  value  of  B  is  determined 
and  substituted  above  to  yield  a  Larmor  radius  of  3.7 MO'4  [m].  This  is  ~  one  fourth  of  the 
assumed  pinch  radius  of  0.0015  [m]. 
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From  a  qualitative  standpoint,  then,  the  diamagnetic  effect  is  quite  dominant,  as  the  Lai  mor 
radii  of  the  pinch-edge  deuterons  show  that  they  gyrate  into  the  pinch  considerably.  Away  from 
the  pinch  center,  toward  the  ends  of  the  pinch,  the  effect  is  more  noticeable  because  the  Larmor 
radii  comprise  a  larger  radial  fraction  of  the  pinch  volume. 

In  developing  a  two-dimensional  code,  a  traditional  cylindrical  coordinate  system  was 
chosen,  for  the  pinch  exhibits  approximately  cylindrical  shape.  For  the  shape  chosen  to  model  the 
pinch  with,  the  only  digression  from  a  cylindrical  shape  is  the  narrowing  of  the  pinch  at  the  axial 
ends  of  the  pinch  region,  thus  resulting  in  the  parabolic  contour  at  the  pinch  edge.  For  one  half  of 
the  pinch,  examples  of  the  assumed  "streamlines"  of  particles  flowing  in  the  pinch  are  shown  in 
Figure  3.3.  The  streamlines  should  be  considered  as  computational  grid  tools  only,  for  they  are 
used  to  represent  the  general  velocity  of  a  group  of  particles  around  a  particular  grid  point.  The 
"streamlines"  do  not,  then,  necessarily  represent  the  actual  paths  of  particles  through  the  pinch,  as 
the  pinch  is  assumed  to  be  in  equilibrium.  In  this  case  it  is  assumed  that  no  particles  leave  the 
pinch,  save  those  that  are  products  of  thermonuclear  fusion  within  the  pinch.  To  determine  the 
equilibrium  profiles  using  the  numerical  grid  computation,  it  is  necessary  to  use  the  correct  value  of 
the  azimuthal  magnetic  field  induced  by  the  plasma  current  at  each  axial  division  of  the  pinch 
region. 

Here  is  where  the  results  of  diamagnetism  come  into  play.  It  is  of  key  importance  to 
understand  to  what  degree  the  azimuthal  magnetic  field  is  assumed  nullified  by  the  diamagnetic 
field  induced  by  the  magnetic  current  flowing  in  the  pinch.  If  its  effect  is  near-complete  dominance 
within  a  skin  layer  at  the  edge,  then  the  magnetic  field  within  the  pinch  volume  is  considered  to  be 
near  zero.  For  the  computations  used  to  produce  the  results  in  this  report,  the  azimuthal  magnetic 
field  is  assumed  to  fall  off  in  linear  fashion,  reaching  zero  at  the  pinch  axis.  An  illustration  of  this 
magnetic  field  degradation  is  shown  in  Figure  3.4,  where  the  Be  term  is  being  driven  to  zero  by  the 
diamagnetic  current  in  the  plasma.  Where  the  illustration  shows  the  field  to  degrade  sharply  at 
first,  then  more  gradually  approaching  the  zero  value,  the  assumption  here  is  that  the  degradation 
exhibits  a  straight-line  behavior.  A  20  by  20  grid  was  used  for  the  pinch  region  model. 

Fusion  power  density  calculations  were  done  using  the  averaged  values  of  ion  number 
densities  and  using  a  reaction  rate  parameter  value  based  on  the  average  kinetic  temperature  within 
the  pinch.  The  fusion  power  density  is  given  by: 


PDfusion  =  n,n2«Tv>Ei2  ,  (3.10) 

where  nj  and  n2  are  the  number  densities  of  the  two  reactant  species,  <av>  is  a  term  known  as  the 
reaction  rate  parameter,  or  average  product  of  the  thermonuclear  cross  section  a  and  the  relative 
speed  v  between  species  1  and  2,  and  E12  is  the  energy  yield  per  thermonuclear  reaction  of  one  pair 
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Figure  3.3 

Streamline  Map  for  Half  of  the  Pinch  Region  for  Radial  vs.  Axial  Distances 
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Figure  3.4 

Schematic  of  Diamagnetic  Degradation  of  Azimuthal  B  Field  [31] 
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of  reactant  species.  Values  of  the  reaction  rate  parameter  are  well-documented  in  various  texts 
[25],  as  are  the  energy  yields  per  reaction,  obtained  through  a  simple  conservation  of  mass/energy 
calculation  for  the  appropriate  thermonuclear  reaction.  The  energy  yields  for  the  D-D  and  D- He 
reactions  are  3.65  [MeV]  and  18.3  [MeV],  respectively,  where  the  value  for  D-D  is  the  average 
energy  yield  of  the  two  major  reaction  pathways,  which  have  approximately  equal  occurrence 
probabilities.  Since  <ov>  is  a  function  of  the  energy  kT  of  the  pinch  particles,  the  profiles  of  n 
and  kT  obtained  in  the  next  section  on  pinch  scaling  relations  can  be  used  with  available  <cv>  data 
and  Eq.  (3.10)  to  yield  the  fusion  energy  density  output  of  the  pinch 

Assuming  the  pinch  to  have  a  circular  cross  section,  which  follows  from  assuming 
azimuthal  symmetry,  the  volume  of  the  pinch  is  equivalent  to  V  =  Jt(rav2)(/pinch)-  Thus,  a  typical 
pinch  volume  would  be:  1.50-1  O'7  [m3].  If  a  D-3He  fuel  mixture  is  considered,  with  np  =  nHi;-3  = 
0.5n,  then  Eq.  (3. 10)  becomes: 

1  ,  1  2 

PDf  =  — <n>  <(<ov>D  3Hc)>^D3He  +  g’<n>  <(<<Jv>Dd)>EdD  ’  (3.11) 

where  the  second  set  of  averaging  signs  surrounding  each  <ov>  and  each  n  value  indicates  an 
average  over  the  pinch  region,  both  radially  and  axially.  Such  averages  are  computed  during  the 
execution  of  the  numerical  computation,  just  as  the  <r>  value  discussed  above  is  computed.  For 
energies  in  the  range  of  -  i  to  100  [keV],  the  <crv>  values  are  given  by  the  log  base- 10  curve  fits 
found  in  [25].  The  kT  values  used  in  the  equations  found  there  are  input  in  units  of  [keV],  while 
the  resulting  values  of  <av>  are  in  units  of  [cm3/s].  Keeping  with  the  MKS  system  used  in  this 
report,  the  <av>  values  are  multiplied  by  a  factor  of  106  to  obtain  units  of  [m3/s].  The  fusion 
power  output  is  computed  by  multiplying  the  fusion  power  density  by  the  volume  of  the  pinch. 
From  this,  a  graphical  representation  can  be  made  of  the  fusion  power  versus  pinch  current  to 
determine  a  scaling  relation. 
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PINCH  SCALING  RELATIONS  AND  THEIR  APPLICATION  TO  THE  DPF  PROPULSION 
SYSTEM 


Review  of  Pinch  Scaling  Relations 

Experiments  at  the  University  of  Illinois  have  investigated  the  scaling  of  pinch  dimensions 
in  a  dense  plasma  focus  (DPF)  device  [30].  For  a  deuterium  filling  pressure  of  3.0  [torr],  the 
current  in  the  plasma  was  found  to  have  little,  if  any,  effect  on  the  pinch  length  (Figure  4.1). 
However,  with  regard  to  the  radius  of  the  pinch,  a  small  dependence  was  observed  (Figure  4.2). 
Concerning  the  scaling  of  the  pinch  radius  and  length  with  the  filling  pressure,  the  length  scaling 
was  not  obvious  due  to  the  large  spread  in  data  (Figure  4.3),  while  the  radius  scaling  was  quite 
well-behaved  (Figure  4.4),  as  compared  to  the  three  afore-mentioned  scaling  plots.  For  each  of  the 
four  observed  relationships,  the  power  fit  law  is  stated  on  the  appropriate  figure  from  Figures  4.1- 
4.4.  In  each  case  the  exponent  value  is  considerably  less  than  unity;  perhaps  even  exhibiting  zero 
dependence  on  the  property  in  question.  These  experimental  results  show  the  most  recent  work 
found  regarding  experimentally  determined  pinch  size  scaling  studies  of  the  DPF. 

In  the  research  described  in  this  report,  differing  values  of  pinch  radius  and  length  were 
specified  to  obtain  the  temperature  and  number  density  profiles,  and  currents  associated  with  each 
set  of  values.  It  was  also  possible  to  change  the  given  values  of  maximum  number  density  and 
drift  velocity  in  order  to  obtain  how  the  fusion  power  output  scaled  with  the  plasma  current.  The 
number  density  profile  for  tokamak-like  plasmas  was  assumed  in  the  radial  direction  of  the  pinch, 
with  the  maximum  value  of  number  density  occurring  at  the  center  of  the  pinch,  while  along  the 
axial  direction  of  the  pinch,  a  modified  form  of  the  tokamak  profile  exhibiting  a  lesser  gradient  was 
used  in  order  to  agree  with  the  axial  profiles  given  by  Potter  [32]. 

For  comparison,  use  was  made  of  work  done  by  Potter  [32]  in  order  to  see  how  the 
average  ion  number  density  (n),  average  ion  kinetic  temperature  (kT),  average  fusion  power 
density,  and  average  axial  velocity  (vz)  vary  over  the  length  of  the  pinch  region,  per  Figure  4.5 
showing  typical  results  from  the  present  work.  Using  a  two-dimentional  2()-by-2()  mesh  (Figure 
4.6)  and  the  plasma  ion  axial  current  based  on  the  given  ion  number  density  profile  (Figure  4.7) 
and  the  axial  velocity,  the  local  magnetic  field  profile  (Figure  4.8)  and  the  kinetic  temperature 
profile  (Figure  4.9)  were  determined.  The  B  profile  follows  the  earlier  conservation  relation  set 
forth  in  Eq.  (3.6).  From  the  appearance  of  the  kT  profile,  it  can  be  seen  that  the  radial  profile  has 
the  same  general  behavior  of  the  n  profile  which  was  assumed  to  be  tokamak-like.  The  only 
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Figure  4.2 

Plasma  Current  Scaling  of  Pinch  Radius  [30} 
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Deuterium  Filling  Pressure  Scaling  of  Pinch  Radius  [30] 
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Figure  4.5  _ 

Typical  Test  Axially  Averaged  Profiles  of  Ion  Number  Density  n.  Kinetic  Temperature  kT,  Axial 

Velocity  v_z,  and  Fusion  Power  Density 
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Figure  4.6 

Computational  Mesh  Used  in  Pinch  Calculations  (20-by-20) 
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digression  from  the  tokamak  shape  for  the  kT  profile  is  that  the  profile  approaches  maximum 
values  much  more  rapidly  moving  into  the  pinch  radially  than  does  the  n  profile.  Thus,  the  factor 
in  the  kT  profile  may  be  more  appropriately  given  by  0.1  rather  than  the  assumed  0.9  of  the 
tokamak  profile  used  in  establishing  the  m  =  0  drift  model's  feasibility  (see  Eq.  (3.1)).  This  is  also 
evidenced  by  the  kT  contour  lines  (solid  lines)  appearing  in  Figure  4.9,  while  the  streamlines  in  the 
pinch  are  marked  by  the  dashed  lines,  which  follow  the  shape  of  the  pinch.  By  varying  the 
current,  a  number  of  runs  yielded  the  corresponding  thermonuclear  reaction  power  yields,  giving 
the  relation  illustrated  in  Figure  4.10,  from  which  a  scaling  of  Pf  ~  Ip5  4  is  found.  Figure  4.1 1 
shows  that  the  inherent  current-squared  scaling  of  kinetic  temperature  density  (or  energy  density) 
as  assumed  in  this  study  holds.  Closer  investigation  will  reveal  a  slope  of  2  on  the  log-log  plot  of 
Figure  4.1 1.  One  should  bear  in  mind  that  Figures  4.10  and  4. 1 1  are  log-log  plots.  Calculating 
the  slopes  on  such  plots  yields  the  scaling  exponent  directly. 

Now  earlier  it  was  established  that  experimental  scalings  of  neutron  yield  with  current  were 
characterized  by  an  exponent  value  of  3.3  to  4.0  for  the  pinch  current  in  the  scaling  relation.  Of 
course,  such  relations  were  based  on  the  use  of  a  D-D  fuel  mixture,  while  here  we  are  utilizing  a  D- 
3He  fuel  mixture.  The  use  of  this  fuel  mixture  in  experiments  has  been  rare,  but  a  group  led  by 
Hirano  did  make  use  of  it  [33].  They  used  the  ratio  of  D-D  neutrons  produced  to  D-3He  protons 
produced  to  determine  the  temperature  of  the  plasma.  In  the  kinetic  temperature  (kT)  region  of 
interest  for  the  runs  which  were  performed  (typically  average  kT  values  between  8  and  15  keV), 
the  reaction  rate  parameter  curve  for  D-3He  varies  more  rapidly  than  the  D-D  curve,  thus  offering  a 
reason  for  the  scaling  exponent  of  the  pinch  current  being  larger  than  those  reported  for  D-D 
systems.  Because  fusion  power  is  directly  proportional  to  the  reaction  rate,  the  fusion  power 
should  also  be  directly  proportional  to  the  number  of  particles  produced  for  a  given  reaction. 
Therefore,  it  is  the  reaction  rate  parameter  curves  which  are  responsible  for  the  favorable  scaling  of 
the  fusion  power  with  pinch  current 


DPF  Propulsion  System  Application  of  Scaling  Relations 

From  the  findings  discussed  in  the  previous  section  on  pinch  scaling  relations,  the 
application  of  such  results  was  introduced  into  a  DPF  propulsion  code  developed  earlier  [1,2,23]. 
Before  relating  the  results,  a  brief  explanation  of  appropriate  performance  characteristics  for  a 
space  propulsion  system  is  warranted.  Continuing  from  a  discussion  of  fusion  power  output 
calculations  mentioned  previously  in  pinch  scaling  relations,  one  remembers  that  the  fusion  power 
output  is  equal  to  the  product  of  fusion  power  density  and  pinch  volume.  The  usable  fraction  of 
this  power  is  reduced  by  radiation  loss  power  in  the  form  of  bremsstrahlung  and  synchrotron 
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radiation  [26],  The  remaining  power  recharges  the  capacitor  banks  on  board  the  DPF  propulsion 
system  and  heats  the  propellant  released  during  its  operation. 

The  portion  which  heats  the  propellant  is  the  exhaust  power  Pcx.  The  combination  of 
reaction  products  and  propellant  spewing  forth  from  the  exit  of  the  DPF  system  creates  the  thrust 
which  propels  the  vehicle.  The  exhaust  velocity  of  the  propellant  is  written  in  terms  of  the  exhaust 
power  and  the  exhaust  thrust  Fex  as  [23]: 


v 


ex 


(5.1) 


where  the  factor  of  2.0  has  come  from  exhaust  velocity  enhancement  due  to  the  use  a  magnetic 
nozzle  to  direct  the  exiting  fusion  reaction  products  and  ionized  propellant  particles.  From  the 
value  of  vex,  one  of  the  primary  performance  parameters  is  determined-the  specific  impulse  (Isp). 
Roughly,  Isp  can  be  considered  as  the  efficiency  of  the  vehicle's  propulsion  system  as  a  whole, 
i.e.,  its  "miles  per  gallon"  value.  High  values  of  Isp  can  enable  missions  requiring  interplanetary 
traveling  distances  and  large  changes  in  vehicle  velocity  that  systems  with  low  specific  impulses 
cannot  achieve.  For  a  thrust  parallel  to  the  exhaust  velocity  (Fex-vex  =  Fexvex),  the  specific 
impulse  of  the  system  is  written  as: 
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sp 


(5.2) 


and  the  thrust  by: 


^propeilantvex  ■  (5.3) 

Here,  mpropeiiant :  j  the  mass  flow  rate  leaving  the  vehicle.  Of  the  mass  flow,  the  fraction  of  it 
which  is  propellant  is  much,  much  greater  than  the  portion  of  it  which  is  fusion  product  particles, 
so  the  mass  flow  rate  is  approximated  by  the  mass  flow  rate  of  the  propellant  alone  with  very  little 
error.  The  term  g  is  simply  the  gravitational  acceleration  of  the  earth  at  sea  level.  With  g  as  the 
chosen  constant  of  acceleration,  Isp  is  used  as-  a  gauge  of  performance  to  compare  space  propulsion 
systems. 

Rather  than  the  thrust  alone  being  used  as  a  benchmark  for  space  vehicle  performance,  a 
quantity  of  greater  importance  is  the  thrust-to-weight  [F/W]  ratio,  which  is  simply  the  ratio  of  the 
thrust  F  to  the  total  weight  of  the  vehicle  (W  =  mveh  g)  at  a  given  time,  i.e.,  the  system  "horse 
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power."  Use  of  the  thrust-to- weight  ratio  and  the  specific  impulse  together  provides  a  good  way 
of  comparing  space  propulsion  system  performances. 

Utilizing  the  resulting  profiles  obtained  in  this  study,  the  existent  propulsion  system  code 
was  altered  to  include  the  detailed  pinch  calculations  considered  in  previous  sections.  For  a  one- 
thruster  base  case  of  Av  =  10  [km/s],  pinch  current  =  20  [MA],  and  propellant  mass  flow  rate  =  30 
[kg/s],  the  major  performance  parameters  listed  in  Table  IV  are  obtained.  Also  shown  are  the 
results  given  in  the  earlier  study  of  enhanced  electrode  designs  [2].  In  that  work,  the  parameters 
shown  in  Table  IV  were  obtained  for  this  base  case  considering  the  use  of  the  Livermore-I  DPF 
and  then  the  use  of  the  enhanced  electrode  design  (EED)  DPF  which  was  determined.  It  is  evident 
that  a  more  detailed  treatment  of  the  pinch  region  has  led  to  less  promising  values  of  the  F/W  ratio 
and  specific  impulse  for  the  given  pinch  current  of  20  [MA],  although  at  higher  pinch  currents  of 
30  to  40  [MA],  the  more  detailed  treatment  produces  results  perhaps  even  more  favorable  than 
those  mentioned  in  earlier  reports. 


Table  IV 

Comparison  of  Performance  Parameters  for  Livermore-I  DPF, 

Enhanced  Electrode  Design  DPF,  and  Detailed  Pinch  DPF  Propulsion  System  Test  Cases 


Parameter  Livermore-I  DPF 

Total  Mass  [kg]  3.69E5 

Total  Thrust  [N]  4.73E5 

Thrust-to-Weight  0. 1 3 1 

Specific  Impulse  [s]  1583 


EED  DPF 
3.59E5 
4.84E5 
0.137 
1607 


Detailed  Pinch  DPF 
5.59E5 
9.07E4 
0.0165 
973.6 


After  running  various  cases  with  the  new  scaling  considerations  and  pinch  calculations 
incorporated  into  the  existing  propulsion  code  [1],  some  typical  mission-required  axial  profiles  are 
shown  in  Figure  4. 12.  In  comparing  these  profiles  with  the  earlier  test  case  profiles  of  Figure  4.5, 
it  is  evident  that  the  parameters  required  in  space  propulsion  are  greater  value.  In  the  two  cases 
being  compared,  the  axial  velocity  distribution  is  the  same,  while  the  number  density  profile  has 
been  increased  by  an  order  of  magnitude.  The  resulting  kT  axial  profile  shows  the  same  behavior, 
but  its  value  has  increased  an  order  of  magnitude  as  well.  What  is  striking  is  the  considerable 
increase  in  the  fusion  power  density  emerging  from  the  pinch-an  increase  of  over  4  orders  of 
magnitude.  The  distribution  of  power  emerging  from  the  pinch  has  also  become  less  centralized 
axially.  For  illustration,  the  assumed  ion  number  density  profile  for  a  typical  mission-required 
regime  of  operation  is  shown  in  Figure  4.13,  while  the  subsequently  determined  profiles  of  local 
magnetic  field  and  kinetic  temperature  are  presented  in  Figures  4.14  and  4.15,  respectively. 

In  order  to  gain  an  understanding  of  the  effects  on  propulsion  performance,  a  set  of  18  test 
cases  were  considered,  in  which  the  basic  variables  were  the  axial  velocity  in  the  pinch  and  the  ion 
number  density  in  the  pinch.  Figure  4.16  gives  a  surface  representation  to  each  of  these  18  test 
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Figure  4.12 

Typical  Mission-Required  Axially  Averaged  Profiles  of  Ion  Number  Density  n.  Kinetic 
Temperautre  kT,  Axial  Velocity  v_z,  and  Fusion  Power  Density 
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Figure  4.13 

Typical  Mission-Required  Input  Ion  Number  Density  Profile 
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Figure  4.16 

Pinch  Currents  Associated  with  the  18  Test  Case  Combinations  of  Ion 
Number  Density  and  Axial  Velocity 


cases.  Three  values  of  ion  number  density  were  used:  1.0,  2.0,  and  3.0- 1026  [m*3],  while  six 
values  for  the  pinch  axial  velocity  were  chosen:  0.7,  0.9,  1.0,  2.0,  3.0,  and  4.0- 10s  [m/s].  In 
choosing  these  values,  the  previous  project  work  regarding  the  propulsion  code  was  considered,  as 
well  as  observances  of  the  results  obtained  from  making  a  number  of  test  runs  before  deciding  on  a 
precise  test  grouping.  By  observing  the  I  =  0  plane  of  Figure  4. 16,  the  projections  of  the  test  case 
points  from  the  surface  are  visible.  It  was  generally  noted  that  velocities  below  0.7- 105  [m/s] 
could  not  generate  enough  fusion  power  to  heat  any  propellant  mass,  while  number  densities  and 
axial  velocities  beyond  those  chosen  yielded  unheard  of  values  for  the  pinch  current,  meaning 
currents  on  the  order  of  100  [MA]  or  greater. 

Some  general  findings  will  now  be  presented.  Figure  4.17  presents  a  base  set  of 
parameters  including  a  4-thruster  device  and  a  mission  velocity  increment  of  10  [km/s].  By 
varying  the  propellant  mass  flow  rate  from  10  to  20  to  30  [kg/s],  it  was  shown  that  increasing  the 
rate  increases  the  F/W  ratio  and  decreases  the  specific  impulse.  This  remains  consistent  with  the 
results  found  in  earlier  project  reports  [1,2],  as  expected.  The  vast  majority  of  the  output  found  in 
this  report  is  akin  to  that  shown  in  Figures  4.18  and  4.19,  where  for  specified  values  of  n,  vz,  Av, 
and  propellant  mass  flow  rate,  the  resulting  values  of  pinch  current,  specific  impulse,  and  thrust- 
to-weight  ratio  are  plotted  as  a  surface  for  a  number  of  thrusters  varying  from  1  to  4.  In  total,  five 
values  of  propellant  mass  flow  rate  (1,  5,  10,  20,  and  30  [kg/s])  and  three  values  of  Av  (10,  20, 
and  40  [km/s])  were  chosen  to  vary  with  the  1,2,  3,  and  4-thruster  configurations.  This  results  in 
60  total  mission  cases,  not  to  be  confused  with  the  18  test  cases  previously  mentioned.  Only  the  1 
[kg/s]  flow,  10  [km/s]  Av  cases  are  shown  in  the  present  section,  with  the  remainder  of  them 
found  in  Appendix  A. 

In  Figure  4.20  a  surface  plot  of  the  variance  of  DPF  system  operating  energy  as  a  function 
of  the  F/W  ratio  and  the  specific  impulse  is  presented  for  the  purpose  of  visualizing  the  energy 
requirements  of  such  a  system,  while  in  Figure  4.21,  scaling  laws  for  the  fusion  yield  are 
revisited.  It  was  observed  how  the  key  performance  parameters  for  the  18  test  cases  varied  when 
considering  fusion  number  yield  scaling  of  pinch  current  to  the  2nd,  3rd,  4th,  and  5th  powers.  It 
is  evident  that  for  a  given  test  case,  the  lower  scaling  is  more  favorable,  at  least  for  the  higher 
performance  cases,  coinciding  with  the  30  to  40  [MA]  "magic"  zone  of  operation  before 
performance  begins  falling.  This  magic  zone  behavior  is  present  in  all  of  the  plots  of  the  type 
found  in  Appendix  A.  Physically,  it  seems  that  the  culprit  behind  the  lessened  performance  is  the 
reaction  rate  parameter  curve  for  the  D-3He  reactions,  which  begins  to  plateau  and  decrease  for 
kinetic  temperatures  exceeding  100  [keV].  Also,  regarding  the  seemingly  better  performance  for 
lower  scalings,  it  should  be  noted  that  a  fuel  cycle  modeling  is  not  present  in  the  propulsion  code. 
Once  a  fusion  fuel  particle  has  reacted,  it  is  considered  to  not  have  any  effect  on  subsequent  fusion 
reactions.  The  decreasing  number  density  and,  thus,  the  decreasing  current  over  the  lifetime  of  the 
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Figure  4.17 

Effect  of  Propellant  Mass  Flow  Rate  on  Specific  Impulse  and  FAV  Ratio  in  18  Test  Cases  at  3 
Different  Rates  for  a  Av  of  10  [km/s]  Using  a  4-Thruster  System 


(a)  1  Thruster  (+’s  in  F/W=0  plane) 
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Figure  4.18 

Pinch  Current,  Specific  Impulse,  and  F/W  Ratio  for  Propellant  Mass  Flow  Rate  of  1.0  [kg/s]  and 
Av  =  10  [km/s]  with  (a)  l  Thruster  and  (b)  2  Thrusters 
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Operating  Energies  for  Test  Cases  at  3  Propellant  Flow  Rates  Using  4  Thrusters 
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Figure  4.20 

DPF  Operating  Energies  vs.  F/W  Ratios  and  Specific  Impulse  Values  for  Propellant  Flow  Rates  of 
10, 20,  and  30  [kg/s]  Using  4  Thrusters  with  a  Av  of  10  [km/s] 
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Figure  4.21 

Effect  on  Performance  Parameters  for  Various  Scalings  of  Fusion  Number  Yield  with  Pinch 

Current  (Yp  ~  I“,  with  n  =  2,3,4,5) 
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pinch  are  more  noticeable  at  the  higher  scaling  exponents.  It  would  seem  that  the  quicker  depletion 
of  the  resource  fuel  particles  early  in  the  pinch  phase,  while  leading  to  higher  initial  power  output, 
also  leave  fewer  particles  at  the  end  of  the  pinch  phase  to  react.  Thus,  the  lifetime  of  the  pinch 
seems  to  also  play  a  key  role  in  the  scaling  considerations  for  the  pinch,  at  least  as  far  as 
propulsion  performance  is  concerned. 
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CONCLUSION  AND  RECOMMENDATIONS 


In  completing  this  current  three-year  project,  work  has  moved  from  an  overall  propulsion 
system  standpoint  inwardly  to  the  mechanisms  which  are  at  work  in  the  pinched  plasma  itself.  It 
was  determined  during  the  first  year's  work  that  an  impulsive  mode  of  operation  gave  the  best 
propulsion  system  performance.  For  a  mission  requiring  a  Av  of  40  [km/s],  a  F/W  ratio  of  almost 
0.08  and  a  specific  impulse  of  4000  Is]  were  possible  with  4  thrusters  operating  at  20  IMA).  This 
was  the  optimum  operating  regime  for  the  impulsive  mode.  Continuous  firing  modes  fell  two 
orders  of  magnitude  short  of  the  required  F/W  ratios  for  manned  space  travel.  Also  researched  in 
the  first  year  of  the  project  was  the  benefits  derived  from  spin  polarization  of  the  D-3He  fusion 
fuel.  Power  output  increases  from  the  fusion  reactions  of  20%  to  50%  were  found  to  be  possible 
in  relation  to  the  reaction  rates  of  the  thermonuclear  reactions.  This  consequently  increases  both 
the  F/W  ratio  and  the  specific  impulse  of  the  system.  Spin  polarization  may  also  work  to  suppress 
the  D-D  side  reactions,  thus  meaning  that  the  power  increase  of  50%  could  be  considered  a  total 
power  increase  in  the  system  [1]. 

Moving  into  the  second  year  of  work  in  the  project,  the  effects  of  the  electrical  system 
configuration  on  the  propulsion  performance  figures  of  merit  were  investigated.  An  equivalent 
circuit  of  the  DPF  device  was  devised  for  use  in  a  1-D  transient  code  to  determine  sheath  dynamics 
during  the  rundown  phase.  The  leakage  current  was  modeled  as  a  constant  resistance  which 
removed  current  from  the  plasma  sheath.  The  Livermore-I  DPF  was  used  for  comparison  with  the 
test  cases  thus  obtained.  Various  radial  and  axial  anode  variations  were  studied  as  a  means  of 
modifying  the  Livermore-I  device  to  achieve  the  necessarily  high  current  needed  to  heat  the  pinch. 
It  was  found  that  the  electrodes  did  not  need  to  undergo  formidable  changes  in  order  to  achieve  the 
high  current  requirements  of  10  to  20  [MA].  Instead,  it  was  the  charging  voltage  which  needed 
drastic  increases,  with  those  increases  being  greater  than  what  is  available  today.  For  the  state  of 
the  art  in  high  current  pulse  power  technology,  the  reference  is  the  SHIVA  experiment  at  Kirtland 
Air  Force  Base,  which  utilizes  a  load  current  of  9  [MA]  at  a  charging  voltage  of  125  [kV],  An 
increase  of  -25%  in  the  charging  voltage  would  be  required,  as  well  as  a  circuit  which  could 
handle  the  high  current  levels.  The  end  result  of  altering  the  electrode  configuration  is  a  lessening 
of  the  required  capacitor  bank  mass,  which  necessarily  increases  the  F/W  ratio.  The  study 
provided  a  more  realistic  view  of  the  requirements  needed  for  space  propulsion  feasibility  of  the 
DPF,  but  the  study  also  provided  a  more  optimistic  view  of  DPF  thruster  performance  as  a  result 
of  the  electrode  modification  investigation  [2], 
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These  first  two  years  of  work  led  into  the  final  phase  of  the  project:  scaling  laws  and 
plasma  pinch  phase  physics.  Such  work  has  been  alluded  to  in  the  conclusions  of  the  past  two 
annual  reports.  Looking  at  the  scaling  of  pinch  size  with  pinch  current  has  shown  that  modeling 
the  pinch  column  as  a  half  wavelength  of  the  m  =  0  instability  based  on  the  electron  drift  velocity 
appears  to  be  feasible,  having  shown  that  number  densities  obtained  in  such  a  model  are  of  the 
correct  order  of  magnitude  for  a  DPF  device.  It  remains  to  be  investigated,  however,  if  this  model 
holds  for  higher  current  configurations,  as  the  model  used  here  is  essentially  indepedent  of  the 
pinch  current.  An  appropriate  dependency  on  the  current  would  need  to  be  determined  in  order  to 
test  the  model's  validity.  From  the  standpoint  of  the  fusion  yield  scaling  with  pinch  current,  a 
scaling  of  Pf  ~  Ip5  4  has  been  found  based  on  assumed  tokamak  number  density  profiles  resulting 
from  the  pinch  size  scaling  model.  Pf  is  the  total  fusion  power  produced  in  the  pinch  region,  with 
D-3He  and  D-D  reactions  considered  in  developing  the  reaction  rate  parameter  [25].  Using  the 
assumption  that  the  magnetic  field  penetrates  into  the  plasma  until  reaching  a  value  of  zero  at  the 
axis  of  the  pinch,  the  widely  accepted  Ip2-scaling  of  the  energy  density  (or  kinetic  temperature 
density)  is  obtained.  This  was  done  by  assuming  the  vacuum  magnetic  boundary  pressure  at  the 
edge  of  the  pinch  is  equal  at  all  times  to  the  sum  of  the  particle  and  local  magnetic  pressures. 

In  looking  to  further  research  in  the  pinched  plasma  region,  another  key  parameter  is  the  fill 
pressure  of  the  fusion  fuels.  Further  work  should  be  done  to  see  what  effects  the  pressure  has  on 
the  size  and  shape  of  the  pinched  region,  as  well  as  its  effect  on  trapped  particle  fraction  in  the 
pinch.  From  such  a  study,  one  could  determine  if  a  low  or  high-pressure  mode  of  operation  is 
desirable  for  DPF  operation  from  a  flight  performance  standpoint.  Within  this  idea  of  pressure 
mode  operation,  it  is  desirable  to  know  if  treatment  of  the  plasma  kinetic  temperature  as 
Maxwellian  is  appropriate.  Hirano,  et  al.  [33]  indicate  a  thermonuclear  treatment  may  not  be  viable 
for  their  low-pressure  mode  operation  of  1.5  [torr]  made  up  of  equal  parts  of  deuterium  and 
helium-3.  They  determine  a  temperature  of  1 1  [keV],  while  neutron  time-of-flight  measurements 
give  a  value  of  3  [keV],  thus  owing  to  the  uncertainty  of  the  appropriate  treatment.  Decker,  et  al. 
make  the  assertion  that  if  kT;  <  3kTe,  then  the  measured  neutron  outputs  cannot  be  explained  by 
thermonuclear  processes  [19],  Also,  advances  in  determining  the  fusion  reaction  mechanisms, 
such  as  the  role  of  filamentary  pinch  structures,  could  provide  enlightenment  as  to  the  desired 
pressure  mode  of  operation  [34],  Other  areas  of  possible  research  include  a  study  of  the  effects  of 
electrode  and  insulator  materials  [35],  the  fraction  of  particles  which  can  be  utilized  as  thrust  using 
the  enhancements  from  a  magnetic  nozzle,  and  a  design  of  the  DPF  propulsion  system  and  vehicle 
as  a  whole,  beyond  the  current  schematic  diagrams. 
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APPENDIX  A 


For  the  18  test  cases  mentioned  earlier  in  this  report,  the  graphical  results  for  the  60 
mission  cases  are  illustrated,  continuing  from  the  earlier  presentation  of  these  results  in  Figures 
4.18  and  4.19.  For  the  figures  shown,  2  mission  cases  are  presented  on  each  figure,  while  all  18 
test  cases  are  shown  on  each  figure.  The  "+"  marks  in  the  "F/W  Ratio  =  0"  plane  indicate  the 
projections  from  the  surfaces  shown  regarding  each  of  the  18  test  cases.  Figures  A.  1  and  A.2  are 
Figures  4.18  and  4.19  repeated. 

In  viewing  these  graphs,  one  should  keep  in  mind  that  for  every  four  graphs  relating  to 
cases  of  1, 2, 3,  and  4  thrusters  being  used,  the  values  of  Av  and  mass  flow  rate  are  held  constant. 
Within  each  group  of  four  graphs,  the  parameter  to  watch  is  the  F/W  ratio,  as  the  pinch  current  and 
specific  impulse  will  not  change  among  the  four.  If  comparisons  are  made  to  observe  how  the 
F/W  ratio  and  Isp  values  vary  for  a  set  number  of  thrusters,  then  one  would  need  to  compare  the 
captions  under  appropriate  figures  in  order  to  choose  the  graphs  with  appropriate  values  of  Av  and 
mass  flow  rate.  For  example,  if  one  wishes  to  see  how  the  F/W  ratio  varies  with  Av  for  a  4- 
thruster  system  with  a  mass  flow  rate  of  10  [kg/s],  then  the  graphs  needed  for  the  comparison  are 
Figures  A.14,  A.  16,  and  A.  18.  It  would  be  evident  that  moving  from  a  Av  of  10  [km/s]  to  one  of 
40  [km/s]  causes  the  F/W  ratio  to  decrease  from  above  0.1  to  below  0.02.  The  pinch  current 
values  for  each  of  the  18  test  cases  in  each  graph  remain  unchanged,  as  pinch  current  is  the 
differentiating  characteristic  of  the  test  cases.  Recall  Figure  4.16  to  see  the  assumed  values  of  ion 
number  density  and  axial  velocity  used  to  determine  each  of  the  18  pinch  currents  used. 
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Figure  A.5 
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APPENDIX  B 


Presented  on  the  following  pages  is  the  pinch  calculation  code  used  to  determine  the  kinetic 
temperature  and  local  magnetic  field  profiles  throughout  the  pinch,  as  well  as  the  axially  averaged 
groupings  of  properties.  Imbedded  within  the  code  is  a  subroutine  used  to  determine  the  velocity 
profile  values  at  each  of  the  mesh  points.  Another  FORTRAN  code  was  used  to  create  the  input 
file  containing  the  number  density  profile  values  at  the  mesh  points.  This  input  file  was  the  only 
input  file  used  with  the  pinch  calculation  code.  The  profile  graphs  were  made  using  output  data 
files  from  the  pinch  calculation  code  and  imported  into  MATLAB  to  produce  their  graphical  form. 
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C- 1  4  -  - 

C  Version  created  on  6  April  1 1  ? 4 

C 

C2345-5 

PROGRAM  DPF_RZ_HALF_R?ROF 
IMPLICIT  NONE 

INTEGER.  I. J, K, L, RDI7, ZDI7, I ERR, CCNV 

PARAMETER  (P.DIV  =  40) 

PARAMETER  (ZD IV  =  20) 

DOUBLE  PRECISION  CONST, TP REV, FACTN, FACTT, CURATZ , FACT 
DOUBLE  PRECISION  ZMAT(RDIV+1,  ZDIV+1)  , NONHOM ( RDIV+1 ) 

DOUBLE  PRECISION  ALPHA ( RDIV+1 ) , BETA ( RDIV+1 ) , DR, TFACT 
DOUBLE  PRECISION  SLOPE ( RDIV+1 , ZDIV+1 ), RSPACE ( RDIV+1 , ZDIV+1 ) 
DOUBLE  PRECISION  N (RDIV+1 , ZDIV+1 ), T {RDIV+1 , ZDIV+I ) 

DOUBLE  PRECISION  VZ (RDIV+1 , ZDIV+1 ) , VR ( RDIV+1 , ZDIV+1 ) 

DOUBLE  PRECISION  DVZR (RDIV+1 , ZDIV+1 ) , DVZZ (RDIV+1 , ZDIV+1 ) 
DOUBLE  PRECISION  DVRR  (RDIV+1 ,  ZDIV+1 )  ,  DVRZ  (RDIV+I ,  ZDIV+1 ) 
DOUBLE  PRECISION  RMID, ZMID, RRATIO, BTHETA (RDIV+1 , ZDIV+1 ) 
DOUBLE  PRECISION  MASSD, KEVCON, E, AVTATZ , AVPDZ , TOTP , BEDGE 
DOUBLE  PRECISION  EPS,  NU ,  MUO  ,  RAL ,  AWZZ ,  AVCUR,  AVCURD,  TCURD 
DOUBLE  PRECISION  NMAX , TMAX, PI , VDRIFT, AVNATZ 
DOUBLE  PRECISION  Z, BLOCAL (RDIV/2+1 , ZDTV+1 ) 

DOUBLE  PRECISION  RWEIGH , AVT , TOTT, DZ , SUMAVT, SUMAVN 
DOUBLE  PRECISION  CURENT,  RADIUS,  TOTN,  AVN,  TOTVZ ,  AWZ,  TPOWD 
DOUBLE  PRECISION  KTLOG, RRDHE3 , RRDDN, RRDDP , POWDEN, RAV 
DOUBLE  PRECISION  AVPOWD, VOLUME, EHE3 , EDDN, EDDP, POWFUS 
DOUBLE  PRECISION  OLD ,  RHSIDE , NPF.EV 

PARAMETER  (E  *  1.60219D-19) 

PARAMETER  (KEVCON  =  1.60219D-I6) 

PARAMETER  (MASSD  =  3 . 343468D-27 ) 

PARAMETER  (PI  =  3 . 141592654D0 ) 

PARAMETER  (MUO  =  4 . DO *PI *1 . D-7 ) 

These  energy  outputs  per  reaction  are  in  [J] 

EHE3  =  2 . 93166D-12 
EDDN  =  5 . 23854D-13 
EDDP  =  6.45606D-13 

READ  *,  RMID 
READ  *,  ZMID 
READ  *,  RRATIO 
READ  *,  NMAX 
READ  *,  TMAX 
READ  *,  VDRIFT 
READ  *,  FACT 

EPS  =  ZMID  /  RMID 
NU  =  RRATIO 
RAL  =  RMID  /  NU 


READ  in  the  solution  matrices  with  B.C.'s  and  initial  values 
CALL  BOUNDS (N.T) 

OPEN  (UNIT=1 , FILE= ' output ' , STATUS* ' UNKNOWN' ) 

OPEN  (UNIT=3 , FILE* ' tzpro ' , STATUS* ' UNKNOWN ' ) 

OPEN  (UNIT=4, FILE='nzpro' , STATUS* ' UNKNOWN' ) 

OPEN  (UNIT* 5, FILE*' trpro' , STATUS* ' UNKNOWN ' ) 

OPEN  (UNIT=9 , FILE* ' nrpro ' , STATUS* ' UNKNOWN ' ) 

OPEN  (UNIT=11, FILE* 'tzbeg' .STATUS* 'UNKNOWN' ) 
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o  n  o  on 


cpem  (umit=i2 . file= •  crteg  ■ ,  status =' unknown'  , 

OPEN  (UNIT=21,FILE='bzloc , STATUS =' UNKNOWN '  ; 

OPEN  (UNIT=22  ,  FILE  =  ' brloc '  ,  STATUS =  '  UNKNOWN  '  . 

OPEN  (UNIT=33 , FILE= ' rdim'  . STATUS^ ■ UNKNOWN ' 

OP  EM  (UNIT=44  ,  FILE=  '  Z  '  ,  STATUS*  '  UNKNOWN '  ) 

OPEN  (UNIT* 4 9 , FILE= ' r ' , STATUS =  ' UNKNOWN' ) 

OP EM  (UNIT* 5 4 , FILE* ' b lccp ' , STATUS = ' UNKNOWN' ! 

OPEM  (UNIT=55, FILE='cp' , STATUS* ' UNKNOWN ' ) 

OP  EM  (UNIT=56  ,  FILE*  '  np '  ,  STATUS* '  UNKNOWN '  ) 

OPEM  ( UNIT=60 , FILE* ‘ avt ' , STATUS =' UNKNOWN ' ) 

OPEN  (UNIT=6 1 ,  FILE* '  avn'  ,  STATUS* ' UNKNOWN'  ) 

OPEN  (UNIT=62  ,  FILE=  '  awz  '  ,  STATUS  ='  UNKNOWN'  ) 

OPEN  (UNIT=63 , FILE= ' avpd' , STATUS* ' UNKNOWN ' ) 

OPEN  (UNIT=64 ,FILE='zlen' , STATUS =' UNKNOWN ' ) 

OPEN  (UNIT=70 , FILE* ' zvals ' , STATUS® ' UNKNOWN ' ) 

DO  120  J  =  1 , RDIV/2+1 

DO  110  I  =  1 , ZDIV+1 

z  =  dble(i-l)  *  2.d0*zmid  /  dble(zdiv) 

write  (11,*)  z,  t(j,i> 

write  (12,*)  rspace(j.i),  t(j,i) 

110  CONTINUE 

120  CONTINUE 

C 

C  Fit  the  streamlines  to  the  function  r  =  A  z~2  +  B  z  +  C  and  fill  RSPACE 
C  and  SLOPE  matrices 

CALL  CURVE ( SLOPE , RS  PACE , RMI D , ZMID , RRATIO ) 

DZ  =  2 . DO  *ZMID  /  DBLE(ZDIV) 

CALL  VELOCITY  ( VZ ,  VR ,  DVZR ,  DVZZ ,  DVR R ,  DVRZ ,  SLOPE ,  DZ ,  RSPACE ,  CURENT ,  N , 

&  VDRIFT) 


TCURD  =  O.DO 
RAV  =  0.D0 
SUMAVT  =  0.D0 
SUMAVN  =  0.D0 

DO  500  I  =  1, ZDIV+1 

RADIUS  =  RSPACE (1,1) 

Z  =  DBLE(I-l)  *  2.D0*ZMID  /  DBLE(ZDIV) 

WRITE  (70,*)  Z+0.01D0 
TOTT  =  0.D0 
TOTN  =0.00 
TOTVZ  =  O.DO 
TOTP  =  O.DO 

DO  390  J  =  2, RDIV/2+1 

RWEIGH  =  RSPACE (J, I) 

TOTN  *  TOTN+N ( J, I ) * RWEIGH* (RSPACE (1,1) -RSPACE (2,1) ) 
TOTVZ  =  TOTVZ+VZ(J, I) *RWEIGH* (RSPACE (1, I) -RSPACE (2, I) ) 
390  CONTINUE 

AVNATZ  =  2.DO*TOTN  /  (RSPACE ( 1 , I) ) **2 . DO 
AWZZ  =  2  .  DO  ‘TOTVZ  /  (RSPACE  ( 1 , 1 )  )  **2  .  DO 
C  print  *,  awzz 

CURATZ  =  AVNATZ* AWZZ* E*  PI ‘RSPACE  ( 1, 1)  **2  .  DO 


From  Ampere's  Law. . . (we  get  the  vacuum  8_theta  value) . . . 

BTHETA (1,1)  =  MU 0* CURATZ  /  (2 . DO *PI* RSPACE (1,1) ) 
PRINT  *,  I,  CURATZ,  BTHETA (1,1) 

From  the  beta  =  1  condition  at  the  edge... 
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Since  3_cn.ec a  =  j  sr.  the 
CCMST  -  STHE'T’- 
BTKETA ( RD IV,'  2  *■  1 
T'RDIV  '2+1,  I!  = 
N  ( 1 ,  I )  =  0  .  DO 
T!i. I)  =  0.D0 


-  .  «  »  o 


axes.  .  . 


CONST  .'I .  ?.: 


DO  395  J  =  2.RDIV/2 

3 THETA ( J , I )  =  B THETA i i , I)  - 

&  DBLE ( J- 1 !  /D3LE ! RDIV  2'  *  3THETA -  1 , I 

T  !  J , I )  =  (CONST  -  BTHETA ' J , I i *  *2 . DO '2 . DO / MU  j ;  /  N 
T(J,  I)  =  T ( J,  I)  /  KEVCON 
395  CONTINUE 


DO  400  J  =  2, RDIV/ 2*1 

KTLOG  =  DLOGIO (T(J, I) ) 

RRDHE3  =  0 . 3 53 5 6D0 ’KTLOG’  *3 . DO  -  3 . 3104D0 ’KTLCG* *2 . DO 
&  +  10. 105D0 'KTLOG  -  25.673D0 

RRDHE3  =  l.D-6*10.D0* 'RRDHE3 

RRDDN  =  0 .29811DO*KTLOG**3 . DO  -  2 . 0830D0*KTLOG*KTLCG 
&  *  5 . 70 14D0* KTLOG  -  22.0878D0 

RRDDN  =  1 .D-6* 10. DO** RRDDN 

RRDDP  =  -0. 053051D0*KTLOG**4.D0  +  0 . 57170 *KTLOG* *  3 . DO 
&  -  2 .404 6D0* KTLOG* KTLOG  +  5 . 6436D0 “KTLOG 

&  -  21.971 

RRDDP  =  1 . D-6*10 . DO  ** RRDDP 

POWDEN  =  N(J, I) *N(J,I) * (RRDHE3  *EHE3 /4 . DO  + 

&  RRDDP*EDDN/8 .DO  +  RRDDP  *  EDDP / 8 . DO ) 

R WEIGH  =  RSPACEfJ, I) 

IF  (J.BQ.l)  R WEIGH  =  0.5D0  *  RSPACE(J,I) 

TOTP  =  TOTP+ POWDEN *RWEIGH* (RSPACE (1,1) -RSPACE (2 , I ) ) 
TOTT  =  TOTT+T ( J, I ) “RWEIGH* (RSPACE (1,1) -RSPACE (2,1) ) 

IF  (J.EQ.2)  RAV  =  RAV  +  1 . DO/2 . D0«RSPACE ( 1 , I ) 

TOTN  =  T07N+N ( J, I ) *RWEIGH* (RSPACE (1,1) -RSPACE (2,1) ) 
400  CONTINUE 


AVNATZ  =  2. DO  *  TOTN  /  ( RSPACE ( 1 , I ) ) * *2 . DO 

AVTATZ  =  2. DO  *  TOTT  /  (RSPACE (1 , I ) ) **2 . DO 

SUMAVT  =  SUMAVT  +  AVTATZ ‘AVNATZ 

SUMAVN  =  SUMAVN  +  AVNATZ 

write  (60,*)  AVTATZ 

write  (61,*)  AVNATZ 

write  (62,*)  AWZZ 

AVPDZ  =  2. DO  *  TOTP  /  ( RSPACE (1,1) ) * *2 . DO 
write  (63,*)  AVPDZ 
write  (64,*)  Z 

TCURD  =  TCURD  +  AVNATZ*  AWZZ  *E 
TPOWD  =  TPOWD  +  AVPDZ 
500  CONTINUE 


RAV  =  2. DO  *  DZ  *  RAV  /  <2.D0*ZMID) 
AVCURD  =  2. DO  *  DZ  *  TCURD  /  (2.D0*ZMID) 
AVCUR  =  AVCURD  *  PI  *  RAV  *  RAV 
AVPOWD  =  2. DO  *  DZ  *  TPOWD  /  (2.D0*ZMID) 
AVN  =  SUMAVN  /  DBLE (ZDIV+1) 

AVT  =  SUMAVT  /  DBLE(ZDIV+1)  /  AVN 
VOLUME  =  PI  *  RAV  *  RAV  *  (2.D0*ZMID) 
POWFUS  =  AVPOWD  *  VOLUME 


WRITE ( 1 , * ) 
WRITE ( 1, *) 
WRITE (1, *) 
WRITE ( 1 , * ) 


'Average  Radius  =  ',RAV, '  [m] . ' 

'Average  Power  Density  =  ', AVPOWD/ 1 . D6 , ' 
'Pinch  Volume  =  '.VOLUME,'  [m**3].' 
'Fusion  Power  Produced  =  ', POWFUS/ 1.D6, ' 


[MW/m**3] . ' 
[MW] . ' 
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WRITE !  I ,  ■  ’Average  Pinch  Current  =  ',A7CrJR  1.35 
WRITE  1.*;  'Average  Pinch  Number -Density  =  *  ,  AVM 
WRITE!!,*)  'Average  Pinch  Kinetic  Temperature  = 
WRITE  (1,*)  'Average  Z  Drift  Velocity  =  ' , AWZ ,  ' 

WRITE ( 1 , * i  'Overall  Energy  Density  =  ' , AVN*A7T, ' 

Print  out  the  final  form  of  the  solution  matrices 

DO  600  J  =  1 , RDIV/2+1 

DO  590  I  =  1 , ZDIV+1 

z  =  dble(i-l)  *  2.d0*2mid  /  dble(zdiv) 
write  (4,  *)  2,  n( j, i) 
zmat ( j , i)  =2 

write  (9,*)  rspace(j,i>,  n(j,i) 
write  (3,*)  2,  t(j,i) 
write  (5,*)  rspace(j.i),  t(j,i) 
write  (33,*)  2,  rspace(j.i) 

WRITE  (66,*)  T ( J, I ) 

IF  (J.GT.l)  THEN 

write  (21,*)  z,  blocal(j,i) 
write  (22,*)  rspace(j.i),  blocal(j,i) 

ENDIF 

590  CONTINUE 

write  (44,591)  (zmat ( j , i) , i=l, zdiv+1) 

591  format  (21 ( f 6 . 4 , lx) ) 
write  (44 , * ) 

write  (49,592)  (rspace ( j , i) , i=l, zdiv+1 ) 

592  format  (21 ( f 7 . 5 , lx) ) 
write  (49, *) 

write  (54,595)  (btheta ( j , i) , i=l , zdiv+1 ) 

595  format  (21 ( f 7 . 1 , lx) ) 

write  (54,*) 

write  (55,597)  ( t ( j , i) , i=l , zdiv+1) 

597  format (21 (f6. 2, lx) ) 

write  (55,*) 

write  (56,599)  (n ( j , i) , i=l , zdiv+1) 

599  format (21 (e9 . 3 , lx) ) 
write  (56,*) 

600  CONTINUE 

PRINT  *,  'N  Matrix:' 

DO  700  I  =  1, ZDIV+1 

PRINT  *,  'n  tm-3]  Column  Z  =  ',1 
PRINT*,  (N (L, I) , Lsl , RDIV/2+1 ) 

700  CONTINUE 

PRINT  *,  'T  Matrix: ' 

DO  800  I  a  1, ZDIV+1 

PRINT  *,  'JcT  [keV]  Column  Z  =  '.I 
PRINT*,  (T(L, I) ,L*1, RDIV/2+1) 


800  CONTINUE 

CLOSE 

(UNITal) 

CLOSE 

(UNIT=3 ) 

CLOSE 

(UNIT=4 ) 

CLOSE 

(UNIT=5) 

CLOSE 

(UNIT=9 ) 

CLOSE 

(UNIT=11) 

CLOSE 

(UNIT=12 ! 

CLOSE 

(UNIT=21) 

CLOSE 

(UNIT=22) 

CLOSE 

(UNIT=33 ) 

CLOSE 

(UNIT=44 ) 

CLOSE 

(UNIT=49> 

CLOSE 

(UNIT=54 ) 

,  A  .  *  ,  .  .<6 

r  j  +*%  «  •  1  1  ' 
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on  o  o  o  o  o 


CLOSE  ;UNIT=55) 

close  (ui::t=5o> 
CLOSE  (UNIT=60; 
CLOSE  (UNITa61! 
CLOSE  ( UNIT=52 ) 
CLOSE  (UNIT=63 ) 
CLOSE  (UNIT=64 ) 
CLOSE  ( UNIT=7 0 ) 

STOP 

END 


C23456 

SUBROUTINE  BOUNDS (N,T) 

INTEGER  RDIV, ZDIV 

PARAMETER  (RDIV  =  40) 

PARAMETER  (ZDIV  =  20) 

DOUBLE  PRECISION  N ( RDIV+1 , ZDIV+1 ) , T (RDIV+1 , ZDIV+1 ) 

READ  *,  N 
READ  *,  T 

RETURN 

END 


C23456 

SUBROUTINE  CURVE ( SLOPE, RSPACE, RMID, ZMID, RRATIO) 

IMPLICIT  NONE 

INTEGER  RDIV, ZDIV, I, J 

PARAMETER  (RDIV  =  40) 

PARAMETER  (ZDIV  =  20) 

DOUBLE  PRECISION  SLOPE (RDIV+1 , ZDIV+1) , RSPACE (RDIV+1 , ZDIV+1 ) 
DOUBLE  PRECISION  RMID , ZMID, NU, EPS , Z , RVAL, RAL 
DOUBLE  PRECISION  NUOLD,  EPSOLD,  SITRFR  (ZDIV+1) 

EPS  a  ZMID  /  RMID 
EPSOLD  =  EPS 
NU  =  RRATIO 
NUOLD  *  NU 
RAL  =  RMID  /  NU 

Compute  the  curve  fits  and  load  the  RSPACE  and  SLOPE  matrices 

Load  the  LEFT  END,  and,  therefore,  the  RIGHT  Q4D  also: 

SURFR(l)  =  RAL 
SURFR(ZriVfl)  =  RAL 

DO  1000  J  =  1 , RDIV/ 2 

RSPACE (J,l)  =  RAL  *  DBLE (RDIV/2+1-J)  /  DBLE (RDIV/ 2 ) 
RSPACE (RDIV+2-J, 1)  a  RSPACE(J,L) 

RSPACE (J,ZDI/+1)  a  RSPACE ( J, 1) 

RSPACE (RDIV+2-J, ZDIV+1)  =  RSPACE (J,l) 

1000  CONTINUE 


Load  the  r  =  0  R-AXIS-. 
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ono  onn  non 


1223 


DO  1320  I  =  1 , ZDI7-1 

RSPACE  iRDIY .'2  +  1. I 
SL0?EiRBI7/2+l.I) 
CONTINUE 


=  o.d: 
=  o.  do 


Load  the  z  =  ZMID  Z-A.XIS  ISLOPE  values  only)  : 

DO  1030  J  =  1 ,  P.DI7+1 

SLOPE (J,ZDIV/2+l)  =  0 . DO 
1030  CONTINUE 


Load  the  TOP  and  BOTTOM  boundaries: 

DO  1050  I  =  l.ZDIV-i 

Z  =  2. DO  *  ZMID  •  DBLE(I)  /  DBLE(ZDIV) 

RSPACEU.I  +  l)  =  RVAL(NU,EPS,RAL,Z) 

SURFRt 1+1)  =  RSPACE ( 1 , 1  +  1) 

RSPACE<RDIV+1, 1*1)  =  RSPACE ( 1 , 1+1) 

1050  CONTINUE 

DO  1057  I  =  2 , ZDIV 

SLOPE (1,1)  =  { RSPACE (1,1+1) -RSPACE (1,1-1) )  / 

&  (2. DO  *  2. DO  *  ZMID  *  1.D0  /  DBLE (ZDIV)) 

1057  CONTINUE 


Load  the  remaining  points: 

DO  1080  I  =  2, ZDIV 

DO  1070  J  =  2, RDIV/2 

RSPACE  (J,  I )  =  SURFR(I)  'DBLE (RDIV/2+1-J)  /DBLE  (RDIV/2  i 
RSPACE (RDIV+2-J, I)  =  RSPACE(J.I) 

1070  CONTINUE 

1080  CONTINUE 

DO  1100  I  =  2  ,  ZDIV/2 

DO  1090  J  =  2, RDIV/2 

SLOPE ( J, I )  =  (RSPACE ( J, 1+1) -RSPACE ( J, 1-1 ) )  / 

&  (2. DO  *  2. DO  *  ZMID  *  1.D0  /  DBLE (ZDIV)) 

SLOPE ( J , ZDIV+2 - 1 !  =  -SLOPE (J,  I) 

SLOPE (RDIV+2-J, I)  =  -SLOPE ( J,  I) 

SLOPE (RDIV+2-J, ZDIV+2 -I)  =  SLOPE (J, I ) 

1090  CONTINUE 

1100  CONTINUE 

DO  1107  J  s  2. RDIV/2 

SLOPE (J, 1)  =  (RSPACE ( J ,  2  ) -RSPACE  ( J,  1 )  )  / 

&  (2. DO  *  ZMID  /  DBLE ( ZDIV) ) 

SLOPE (RDIV+2 -J, 1!  =  -SLOPE ( J, 1) 

SLOPE (J, ZDIV+1)  =  -SLOPE ( J, 1 ) 

SLOPE (RDIV+2-J, ZDIV+1)  =  SLOPE (J,l) 

1107  CONTINUE 

SLOPE (1,1)  =  SLOPE (2,1!  +  SLOPE (2 , 1 ) -SLOPE ( 3 , 1 ) 

SLOPE (1, ZDIV+1)  =  -SLOPE (1,1) 

SLOPE (RDIV+1,  1)  =  -SLOPE (1,1) 

SLOPE(RDIV+l, ZDIV+1)  =SLOPE(l,l) 

OPEN  (UNIT=l,FILE='rstuff  ’  ,  STATUS* '  UNXNOWN'  ) 

DO  1141  J  =  1 , RDIV/2+1 
DO  1138  I  =  1 , ZDIV/2+1 

Z  *  2. DO  *  ZMID  *  DBLE(I-l)  /  DBLE (ZDIV) 

WRITE  (1,*)  Z, RSPACE ( J, I) 
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113  3  CCNTIIX’E 

ii-;i  continue 

DO  1147  J  =  1 ,  RDIY,'  2  <•  1 
DO  1145  I  =  l.ZDIV/2+1 

Z  =  2. DO  *  ZMID  *  DBLE(I-l)  /  DELE ( ZDIV; 

WRITE  (1,*)  Z. SLOPS!-*. I) 

1145  CONTINUE 
1147  CONTINUE 

CLOSE  ( UNIT= 1 ! 

RETURN 

END 

C23456 

DOUBLE  PRECISION  FUNCTION  RVAL (NU, EPS , RAL , Z ) 

IMPLICIT  NONE 

DOUBLE  PRECISION  NU.EPS.RAL.Z 

RVAL  =  ( 1 . DO-NU ) / (EPS  *NU )  *  ( Z*Z/ ( EPS*NU*RAL) -2 . DO *Z )  +  RAL 

RETURN 

END 


C23456 

DOUBLE  PRECISION  FUNCTION  DERIV (NU, EPS , RAL, Z) 

IMPLICIT  NONE 

DOUBLE  PRECISION  NU.EPS.RAL.Z 

DERIV  =  2 .DO* (1 .DO-NU) / (EPS*NU)  *  ( Z/ (EPS*NU*RAL)  -  1.D0) 

RETURN 

END 


SUBROUTINE  VELOCITY ( V Z , VR , DVZR , DVZZ , DVRR , DVRZ , 

&  SLOPE, DZ, RSPACE.CURENT.N, VDRIFT) 

IMPLICIT  NONE 

INTEGER  I, J.RDIV, ZDIV 

PARAMETER  (RDIV  =40) 

PARAMETER  (ZDIV  =  20) 

DOUBLE  PRECISION  E, CURENT,VZ(RDIV+1, ZDIV+1) ,VR(RDIV+1, ZDIV+1) 
DOUBLE  PRECISION  PI,DVZR(RDIV+1,ZDIV+I) , DVZZ (RDIV+1 , ZDIV+1) 
DOUBLE  PRECISION  DVRR (RDIV+1 , ZDIV+1) , DVRZ (RDIV+1, ZDIV+1) 

DOUBLE  PRECISION  SLOPE (RDIV+1 , ZDIV+1 ), DZ, RSPACE (RDIV+1, ZDIV+1 ) 
DOUBLE  PRECISION  N(RDIV+1, ZDIV+1) .VDRIFT, FACTOR, VFACT ( ZDIV+ 1 ) 

PARAMETER  (E  =  1.602D-19) 

PARAMETER  (PI  =  3.14159D0) 

C  DATA  VFACT  /0 . 80D0 , 0 . 76D0 , 0 . 71D0 , 0 . 65D0, 0 . 68D0 , 0 . 73D0 , 0 . 77D0 , 

C  &  0 . 83D0 , 0 . 88D0 , 0 . 93D0 , 0 . 97D0 , 1 . 02D0 , 1 . 05D0 , 1 . 07D0 , 

C  &  1 . 08D0 , 1 . 09D0 , 1 . 12D0 , 1 . 15D0 , 1 . 13D0 , 1 . 17D0 , 1 . 20D0/ 

DATA  VFACT  /2 1*1. DO/ 

DO  2000  J  =  1, RDIV+1 

DO  1900  I  =  1, ZDIV+1 
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FACTOR  =  7? ACT i I . 

VR'J.I)  =  SLOPE 'J,  I!  *  '/DRIFT  *  FACTOR 
VZ(J.I)  =  CSQRT(  (VDRI FT "FACTOR/  **2.D0  -  7R  '  J .  I  ;  *  *  2  .  C  : 
C  if i i . eq. 1 ) print  * , r space  ( j , i) , slope ( j , i ) , vz ( j , i ) 

1900  CONTINUE 

2000  CONTINUE 

DO  2200  J  =  2.RDIV 

DO  2100  I  =  2 , ZD IV 

DVRR ( J , I )  =  (VR(J+1, I) -VR1J-1, I) )  / 

&  ( RSPACE ( J+l , I ) -RSPACE ( J- 1 , I ) ) 

DVRZ ( J , I )  =  (VR(J,I+1) -VR(J, 1-1) )  / 

&  (2 . DQ*DZ) 

DVZR(J.I)  =  ( VZ ( J+l , I ) -VZ ( J- 1 , I ) )  / 

4  ( RSPACE ( J+l, I) -RSPACE (J- 1,1) ) 

DVZZ ( J, I )  =  (VZ(J,I+1)-VZ(J,I-1) )  / 

&  ( 2 . DO  *DZ) 

IF  (J. EQ.  RDIV/2  +  1)  THEN 
DVRR (J, I)  =  0.D0 
DVRZ (J, I)  =  0.D0 
DVZR ( J , I )  =  0.D0 
DVZZ (J, I)  =  0.D0 

ENDIF 

2100  CONTINUE 

2200  CONTINUE 

OPEN  (UNIT=2 , FILE= ' velouc ' , STATUS= ' UNKNOWN' ) 

WRITE (2,*)  'DVRR' 

WRITE (2,*)  DVRR 
WRITE (2 , * )  'DVRZ' 

WRITE ( 2 , * )  DVRZ 
WRITE (2,*)  'DVZR' 

WRITE (2,*)  DVZR 
WRITE(2,*)  'DVZZ' 

WRITE (2,*)  DVZZ 
CLOSE  (UNIT=2 ) 

RETURN 

END 
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